# import numpy as np
# from scipy.integrate import solve_ivp
# import matplotlib.pyplot as plt
#
# # 定义非线性多稳态系统的微分方程
# def nonlinear_system(t, state, alpha, beta, epsilon, k):
#     x, y, z, u = state  # 提取状态变量 x, y, z, u
#     # 定义系统方程
#     dxdt = y  # x的导数
#     dydt = x  # y的导数
#     dzdt = -y + 3 * y**2 - x**2 - x * z + alpha + epsilon * u  # z的导数
#     dudt = -k * u - epsilon * (z - beta)  # u的导数
#     return [dxdt, dydt, dzdt, dudt]
#
# # 定义参数
# alpha = -0.02  # 参数 alpha
# beta = 0.15     # 参数 beta
# epsilon = 0.4  # 参数 epsilon
# k = 0.5        # 参数 k
#
# # 初始条件
# initial_conditions = [4.0, 0.0, 0.0, 0.1]  # 初始的 x, y, z, u
#
# # 定义仿真时间
# t_start = 0.0  # 仿真开始时间
# t_end = 60.0   # 仿真结束时间
# t_eval = np.arange(t_start, t_end, 0.01)  # 时间步长为 0.01
#
# # 使用 RK45 方法求解微分方程
# sol = solve_ivp(nonlinear_system, [t_start, t_end], initial_conditions, method='RK45', t_eval=t_eval, args=(alpha, beta, epsilon, k))
#
# # 绘制结果
# plt.figure(figsize=(10, 8))
# plt.plot(sol.t, sol.y[0], label='x(t)')
# plt.plot(sol.t, sol.y[1], label='y(t)')
# plt.plot(sol.t, sol.y[2], label='z(t)')
# plt.plot(sol.t, sol.y[3], label='u(t)')
# plt.title('Nonlinear System Dynamics with RK45')
# plt.xlabel('Time')
# plt.ylabel('States')
# plt.legend()
# plt.grid(True)
# plt.show()
#
# # 将仿真结果保存到npy文件中
# output_data = np.vstack((sol.t, sol.y)).T  # 将时间和状态变量数据组合在一起
#
# # 保存为npy文件
# np.save("nonlinear_system_simulation.npy", output_data)

# import numpy as np
# import matplotlib.pyplot as plt
# from scipy.integrate import RK45
#
#
# # 定义非线性系统的微分方程
# def nonlinear_system(t, y, alpha, beta, epsilon, k):
#     x, y, z, u = y  # 提取状态变量 x, y, z, u
#     dxdt = y
#     dydt = x
#     dzdt = -y + 3 * y ** 2 - x ** 2 - x * z + alpha + epsilon * u
#     dudt = -k * u - epsilon * (z - beta)
#     return [dxdt, dydt, dzdt, dudt]
#
#
# # 使用 RK45 方法求解
# def solve_with_rk45(f, y0, t0, t_end, alpha, beta, epsilon, k, h=0.01):
#     solver = RK45(lambda t, y: f(t, y, alpha, beta, epsilon, k), t0, y0, t_end, max_step=h)
#     t_values = []
#     y_values = []
#
#     while solver.status == 'running':
#         solver.step()
#         t_values.append(solver.t)
#         y_values.append(solver.y)
#
#     return np.array(t_values), np.array(y_values)
#
#
# # 初始条件和参数
# y0 = [4.0, 0.0, 0.0, 0.1]  # 初始条件
# t0, t_end = 0, 60  # 时间范围
# alpha, beta, epsilon, k = -0.02, 0.15, 0.4, 0.5  # 系统参数
# h = 0.01  # 最大时间步长
#
# # 使用 RK45 求解
# t_values, y_values = solve_with_rk45(nonlinear_system, y0, t0, t_end, alpha, beta, epsilon, k, h)
#
# # 绘制结果
# plt.plot(t_values, y_values[:, 0], label='x(t)')
# plt.plot(t_values, y_values[:, 1], label='y(t)')
# plt.plot(t_values, y_values[:, 2], label='z(t)')
# plt.plot(t_values, y_values[:, 3], label='u(t)')
# plt.xlabel('Time')
# plt.ylabel('States')
# plt.legend()
# plt.grid(True)
# plt.title('Solution of Nonlinear System using RK45')
# plt.show()

# import numpy as np
# import matplotlib.pyplot as plt
# from scipy.integrate import RK45
# from matplotlib.animation import FuncAnimation
#
# # 定义非线性系统的微分方程
# def nonlinear_system(t, state, alpha, beta, epsilon, k):
#     try:
#         x, y_var, z, u = state
#         dxdt = y_var
#         dydt = x
#         dzdt = -y_var + 3 * y_var**2 - x**2 - x * z + alpha + epsilon * u
#         dudt = -k * u - epsilon * (z - beta)
#         derivatives = [dxdt, dydt, dzdt, dudt]
#         return derivatives
#     except Exception as e:
#         print(f"Error in nonlinear_system at t={t}: {e}")
#         raise e
#
# # 使用 RK45 方法求解并逐步更新
# def solve_and_plot_realtime(f, y0, t0, t_end, alpha, beta, epsilon, k, h=0.001):
#     solver = RK45(lambda t, y: f(t, y, alpha, beta, epsilon, k),
#                   t0, y0, t_end, first_step=h, max_step=h)
#     print(f"Initial solver status: {solver.status}")
#
#     # 初始化绘图
#     fig, ax = plt.subplots()
#     line_x, = ax.plot([], [], label='x(t)',color='purple')
#     line_y, = ax.plot([], [], label='y(t)')
#     line_z, = ax.plot([], [], label='z(t)')
#     line_u, = ax.plot([], [], label='u(t)')
#
#     plt.legend()
#     plt.grid(True)
#
#     t_values = []
#     x_vals = []
#     y_vals = []
#     z_vals = []
#     u_vals = []
#
#     def init():
#         line_x.set_data([], [])
#         line_y.set_data([], [])
#         line_z.set_data([], [])
#         line_u.set_data([], [])
#         return line_x, line_y, line_z, line_u
#
#     def update(frame):
#         nonlocal solver
#         #print(f"Update called with frame: {frame}, solver status: {solver.status}")
#         if solver.status == 'running':
#             try:
#                 #print(f"Before step: t={solver.t}, status={solver.status}")
#                 solver.step()
#                 #print(f"After step: t={solver.t}, status={solver.status}")
#
#                 t_values.append(solver.t)
#                 x_vals.append(solver.y[0])
#                 y_vals.append(solver.y[1])
#                 z_vals.append(solver.y[2])
#                 u_vals.append(solver.y[3])
#
#                 # 打印当前状态
#                 #print(f"t: {solver.t}, x: {solver.y[0]}, y: {solver.y[1]}, z: {solver.y[2]}, u: {solver.y[3]}")
#
#                 # 更新曲线的数据
#                 line_x.set_data(t_values, x_vals)
#                 line_y.set_data(t_values, y_vals)
#                 line_z.set_data(t_values, z_vals)
#                 line_u.set_data(t_values, u_vals)
#
#                 # 动态调整坐标轴范围
#                 ax.relim()
#                 ax.autoscale_view()
#             except Exception as e:
#                 #print(f"Error during solver step: {e}")
#                 ani.event_source.stop()
#         else:
#             #print(f"Solver status is not 'running' (status: {solver.status}), stopping animation.")
#             ani.event_source.stop()
#         return line_x, line_y, line_z, line_u
#
#     def frame_generator():
#         nonlocal solver  # 添加这一行
#         while solver.status == 'running':
#             #print(f"Frame generator yielding t={solver.t}, solver status: {solver.status}")
#             yield solver.t
#
#     ani = FuncAnimation(fig, update, frames=frame_generator(), init_func=init, blit=False, interval=50)
#     plt.show()
#
# # 初始条件和参数
# y0 = [1.0, 0, 0, 0.1]  # 避免零初始条件
# t0, t_end = 0, 50  # 仿真时间
# alpha, beta, epsilon, k = -0.02, 0.15, 0.4, 0.1  # 调整参数
# h = 0.01  # 更小的时间步长
#
# # 实时绘制解
# solve_and_plot_realtime(nonlinear_system, y0, t0, t_end, alpha, beta, epsilon, k, h)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# 定义非线性系统的微分方程
def nonlinear_system(t, state, alpha, beta, epsilon, k):
    x, y, z, u = state  # 提取状态变量 x, y, z, u
    dxdt = y
    dydt = z
    dzdt = -y + 3 * y ** 2 - x ** 2 - x * z + alpha + epsilon * u
    dudt = -k * u - epsilon * (z - beta)
    return np.array([dxdt, dydt, dzdt, dudt])


# 手动实现 RK45 方法
def rk45_step(f, t, y, h, alpha, beta, epsilon, k):
    k1 = f(t, y, alpha, beta, epsilon, k)
    k2 = f(t + 0.5 * h, y + 0.5 * h * k1, alpha, beta, epsilon, k)
    k3 = f(t + 0.5 * h, y + 0.5 * h * k2, alpha, beta, epsilon, k)
    k4 = f(t + h, y + h * k3, alpha, beta, epsilon, k)
    return y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)


# 使用 RK45 方法求解并逐步更新
def solve_and_plot_realtime(f, y0, t0, t_end, alpha, beta, epsilon, k, h=0.01):
    t = t0
    y = np.array(y0)

    # 初始化绘图
    fig, ax = plt.subplots()
    # ax.set_xlim(t0, t_end)
    # ax.set_ylim(-2, 2)  # 根据情况调整y轴范围
    line_x, = ax.plot([], [], label='x(t)')
    # line_y, = ax.plot([], [], label='y(t)')
    # line_z, = ax.plot([], [], label='z(t)')
    # line_u, = ax.plot([], [], label='u(t)')

    plt.legend()
    plt.grid(True)

    t_values = []
    x_values = []
    y_values = []
    z_values = []
    u_values = []

    def init():
        line_x.set_data([], [])
        # line_y.set_data([], [])
        # line_z.set_data([], [])
        # line_u.set_data([], [])
        #return line_x, line_y, line_z, line_u
        return line_x

    def update(frame):
        nonlocal t, y
        if t < t_end:
            y = rk45_step(f, t, y, h, alpha, beta, epsilon, k)
            t += h

            # 获取当前解
            t_values.append(t)
            x_values.append(y[0])
            y_values.append(y[1])
            z_values.append(y[2])
            u_values.append(y[3])

            # 更新数据
            #line_x.set_data(t_values, x_values)
            # line_y.set_data(t_values, y_values)
            # line_z.set_data(t_values, z_values)
            # line_u.set_data(t_values, u_values)

            line_x.set_data(x_values, y_values)

            ax.relim()
            ax.autoscale_view()

        #return line_x, line_y, line_z, line_u
        return line_x

    def frame_generator():

        nonlocal t  # 添加这一行

        yield t

    ani = FuncAnimation(fig, update, frames=frame_generator, init_func=init, blit=False, interval=50)

    plt.show()


# 初始条件和参数
y0 = [0.2, 0.1, 0.0, 0.1]  # 初始条件
t0, t_end = 0, 50  # 仿真时间
alpha, beta, epsilon, k = -0.02, 0.15, 0, 0.1  # 调整参数，增强非线性影响
h = 0.01  # 最大时间步长

# 实时绘制解
solve_and_plot_realtime(nonlinear_system, y0, t0, t_end, alpha, beta, epsilon, k, h)
